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Nonlinear growth of the bar-mode deformation is studied for a differentially rotating star with 
supercritical rotational energy. In particular, the growth mechanism of some azimuthal modes with 
odd wave numbers is examined by comparing a simplified mathematical model with a realistic 
simulation. Mode coupling to even modes, i.e., the bar mode and higher harmonics, significantly 
f*^ , enhances the amplitudes of odd modes, unless they are exactly zero initially. Therefore, other modes 

^v^j ■ which are not axially symmetric cannot be neglected at late times in the growth of the unstable 

bar-mode even when starting from an almost axially symmetric state. 
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I. INTRODUCTION 

43 : 

There are many three-dimensional simulations that have been carried out for dynamically unstable modes in rotating 
stars for decades QIM,S,IS,[H>[^Q,[S>[S, E3, E3, E3, EE, US, E3, The calculations describe the 

onset of the instability for almost axially symmetric states and its evolution in the nonlinear regime. The critical 
rotational parameter determining the bar- mode instability, i.e., the ratio of the rotational T and gravitational binding 
energies W, has been shown to depend weakly on the rotation law and the equation of state, both in Newtonian 
P, 0, H, 0, HI, @, H H, B EH EH and relativistic gravity [H, El, EH- In addition, recent discoveries from the numerical 
t-H ■ simulations suggest that dynamical bar instabilities can happ en at significantly small values of rotation, which are 
associated with high degree of differential rotation [l5l [Tfjl. IT?! [l8l I20I l2ll I22I |23| . Some of the numerical findings can 
be triggered by the corotation resonance for the so-called low T/W dynamical instability [24], |25[, which is completely 
different in nature from the high T/W one. We here consider the dynamical instability with large rotational parameter, 
high T/W instability only. 

The fate of the unstable bar mode is an interesting problem. The timescale for the persistence of the bar shape is very 
important for the detection of gravitational waves. The problem can only be solved using time-dependent numerical 
codes which require high resolution and accuracy to follow the long-term evolution without being overwhelmed by 
numerical errors. This problem has been independently attacked by various groups using different computational 
codes. Typically, a small non-axisymmetric initial perturbation is added to a rotating axisymmctric model which is in 
unstable equilibrium and the subsequent evolution is calculated numerically. The results show that the bar structure 
is destroyed in a dynamical timescale, i.e., within a few multiples of the rotational period, after the amplitude attains 
;J] ■ nonlinear saturation (see e.g., [TJ EH)- 

In these simulations, some azimuthal modes with odd wave numbers appear in the later growth of the bar shape 
(which has azimuthal number m = 2) (29j . It is, therefore, interesting to study the physical mechanism of the growth 
of odd modes. Are there unstable modes with odd number (e.g., m = 1 or 3) in addition to the unstable bar mode? 
Are the amplitudes of the odd modes enhanced by mode coupling? Three-dimensional simulation of hydrodynamics, 
even in Newtonian gravity, is time-consuming for a wide range of initial data and parameter sets. So the system is too 
complicated for extracting the physical mechanism. In order to understand the growth of odd modes, we investigate 
the evolution of a simplified model. The model's description of mode coupling, unstable growth and decay mimics 
the realistic system very well. Moreover, the number and growth rates of the unstable modes are easily controlled. 
The model, therefore, deepens our understanding of the nonlinear behavior of unstable bar-mode growth in rotating 
stars. The physical mechanism is confirmed by comparing the model problem with a more realistic calculation of a 
dynamically unstable star simulated using three-dimensional hydrodynamics in Newtonian gravity. 
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The organization of this paper is as follows. Recent numerical results for three-dimensional hydrodynamics in 
Newtonian gravity are summarized in Sec. [TTJ A simplified model is presented in Sec. IIIII We investigate how 
odd modes grow within our mathematical model and check that the conclusion is consistent with three-dimensional 
numerical results. Finally a discussion of our results is given in Sec. HVl 



II. THREE-DIMENSIONAL NEWTONIAN HYDRODYNAMICS 

A. Methods 

The set of equations for three-dimensional hydrodynamics in Newtonian gravity is the continuity, Euler, energy and 
Poisson equations: 

d t p + d i (pv i )=0, (1) 

d t (pvi) + d^pvy) = -di(p + p vis ) - pdi® (2) 

dtipe) 1 ^ + ^((p^V) = -I(p £ )-(r-i)/rp vis ^ (3) 

A$ = AnGp. (4) 

Here we assume a T-law equation of state (r = 2), in which thermal pressure P is given by the density p and specific 
internal energy e as 

P = (r - l)pe. (5) 

Pressure P vls in eqs.© and is the artificial viscosity pressure introduced to deal numerically with shocks. The form 
([3]) for the energy equation is derived by eliminating pressure term, and is convenient for the numerical calculation. 
The equations are numerically solved in Cartesian coordinates (x, y, z), assuming planar symmetry across the equator. 
We use PCG(preconditioned conjugate gradient) method e.g, [26| to solve the elliptic equation (g]) with parallel 
processors. The 3D hydrodynamical simulation code in Newtonian gravity has been developed, parallelized and 
tested in the context of dynamical instabilities in Refs. [Til. Il6l. [25j . 

As initial data, we construct differentially rotating equilibrium models with the so-called j-constant rotation law 
with d = 1 (d is a parameter which represents the degree of differential rotation) given by 

ft= m At 2 - (6) 

i?cq + X + V 

Here Q is the angular velocity around z-axis, jo is the constant parameter with units of specific angular momentum, 
and R cq is the stellar radius on equatorial plane. For the construction of the density distribution at equilibrium, we 
also assume a polytropic equation of state with n = 1 

P = K p 2 , (7) 

where n is a constant. 

In the numerical calculations, we monitor the azimuthal Fourier components m — 1, 2, 3, 4 which are defined by 



= T7 [p^^=d s x, (8) 

Q = <e imv ) m=2 

1 f (x 2 - y 2 ) + i[2xy) 3 

= — / p 5 5 a x, (9) 

M J x 2 +y 2 w 

O = (e lmip ) „ 

\ / m— 3 

1 f x(x 2 ~3y 2 ) + iy(3x 2 -y 2 ) ^ 

= Mj P ( X 2 + , 2)3/2 d *> (10) 



M = / im v \ 

' \ / m=4 



1 f (x 4 - 6x 2 y 2 + y 4 ) + i(4xy(x 2 - y 2 )) - 

tt / P ^ « (11) 

M J (x 2 +y 2 ) 2 
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where M is the total rest mass and an angular bracket denotes the density weighted average. 

TABLE I: Four different rotating equilibrium stars in Newtonian gravity of V = 2. 



Model R p /R cq a T/W stability of the bar mode 



I 


0.225 


0.281 


unstable 


II 


0.250 


0.277 


unstable 


III 


0.275 


0.268 


unstable 


IV 


0.300 


0.256 


stable 



a R p : Polar radius; R eq : Equatorial radius 



B. Results for Nonlinear Growth of Non-axially Symmetric Modes 

The growth of the bar-type instability depends on the rotational parameter T/W, the ratio of the rotational 
to gravitational binding energies. The growth rate increases with T/W when it is larger than a critical value of 
T/W « 0.26 ~ 0.27. Our numerical study is limited to the instability with such large rotational parameter. We start 
from hydrostatic equilibrium models for gravity, centrifugal force and pressure gradient. In our numerical experiments 
we found that initial models with very large values of T/W are not suitable for examining the unstable growth, since 
the growing mode too rapidly destroys the bar shape. Models with marginal parameter values are better and we will 
discuss four models (T/W = 0.256 — 0.281) summarized in Table 1. To excite any dynamically unstable mode, we 
disturb the initial equilibrium density p cq with a non-axisymmetric perturbation given by 



£ S^Y m (x,y) 

ro=2,3,4 



(12) 



where S^™ 1 ^ is a small constant, and Y m is a polynomial of order m given by 



y 2 = 4- b 2 + 2x V - V 2 ] > (13) 

Y 3 = -L [x(x 2 - Sy 2 ) + y(3x 2 - y 2 )} , (14) 

^eq 

Ya = ■^ r [x 4 ^6x 2 y 2 +y 4 + Axy(x 2 -y 2 )]. (15) 



cq 



The amplitude of does not coincide with the diagnostics in the azimuthal Fourier components at the initial 
state, rather the relation is given by e.g, |Q| = I^^KJ p C q(x 2 + y 2 )d 3 x) / (V2M R 2 q ) . In order to exclude behaviors 
originating from numerical errors, we check several conditions (the center of mass, conservation of linear momentum 
and conservation of angular momentum) throughout the evolution. All of them are well conserved within several 
percent. We also terminate our code once the relative error in the rest mass exceeds 10~ 4 , because the matter spread 
out from the computational grid. The detailed study for the numerical checks and the accuracy of the numerical 
codes are given in the previous paper 11]. Typically, the integration is terminated at late times of the simulation for 
unstable cases, where turbulent-like behavior is seen. 

The evolution of the bar instability is shown in Figure[T] The initial amplitudes of the perturbation are — 5^ — 
10~ 2 ,<5( 3 ) = 0. It is easily seen that models I-III are unstable, whereas model IV is stable. The exponential growth 
of the m = 2 mode can be seen for models I-III, whereas an oscillation of m = 2 with finite amplitude ~ 5 x 10~ 3 
can be seen in model IV. The m = 4 mode also grows exponentially because of the second harmonic of m = 2 in the 
unstable models (I-III). 

The evolution of the bar instability for different initial perturbation is shown in Figure [5] The rotational parameter 
is the same as that of Figure [TJ but the initial amplitude is reduced to 5^ = 10~ 3 , and that of odd mode is added 
(5( 3 ) = 10~ 3 , = 0. Even for small initial amplitude of m — 2, unstable growth can be seen. The overall features 
of the diagnostics for each component are the same as those of Figure [TJ The unstable bar mode grows in the models 
(I-III), but the saturation time, in which the amplitude of m = 2 mode is ~ 0(1), is delayed due to initially small 
amplitude. In the stable system, model IV, the amplitudes of all diagnostics become same order ~ O(10 -3 ). 
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FIG. 1: Diagnostics as a function of t/P c for four rotating stars I- IV where P c is the central rotation period at t = 0. Solid, 
dashed, dotted, and dash-dotted lines denote the diagnostics for \D\(m — 1), |Q|(m = 2), |0|(m = 3) and |M4|(m = 4), 
respectively. Models I-III are unstable, whereas model IV is stable. 
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FIG. 3: Each four diagnostics as a function of t/P c with four different types of spatial computational resolution. The 
computational grid size is set as (a) (320 x 320 x 80) (solid line), (b) (400 x 400 x 100) (dashed line), (c) (400 x 400 x 100) 
(dotted line), and (d) (500 x 500 x 125) (dash-dotted line), covering the grid points of the equatorial diameter of the equilibrium 
star as (a) 96, (b) 120, (c) 240, and (d) 300. 



The interesting result is that the odd modes grow for late times in the unstable models, irrespective of the initial 
perturbations. The finite differencing scheme used in the numerical code always generates small fluctuations in all m 
modes. This "noise" level is expected to be less than 10~ 4 , or 10~ 3 at most by a rough estimate. That is, the relative, 
numerical error in the global quantities is estimated as 1/N 2 w 10~ 4 -10~ 3 , since second order scheme with typical 
grid number 400 x 400 x 100 is adopted. We also cover the equilibrium star with 120 diameter points in the equatorial 
plane. The amplitudes of odd modes exceed this level in the unstable models (Till) in Figure [TJ The amplitudes of 
m = 3 in Figure [5] exceed it from initial conditions, but they turn to growth after the saturation of unstable mode. 

In order to show that the result is not a numerical artifact with no physics, we have investigated four different types 
of spatial computational resolution. We set the computational grid size as (a) (320 x 320 x 80), (b) (400 x 400 x 100), 
(c) (400 x 400 x 100), and (d) (500 x 500 x 125), covering the equatorial diameter of the equilibrium star with (a) 
96, (b) 120, (c) 240, and (d) 300 grid points. Note that case (a) is 48.8% (= 100 X (1 - 0.8 3 )) reduction of the grid 
points of case (b), and that cases (c) and (d) have twice and 2.5 resolution in the matter regime as that of case (b) 
respectively, keeping the same grid size (the radius of the outer boundary is set closer in cases (c) and (d) than in 
case (b)). We have found that the odd azimuthal modes grow exponentially at the late time in all four cases as shown 
in Figure [3] The initial condition and equilibrium model are the same as those in Figure [TJ-III for cases (a), (b), and 
(c), but we put 1.5 times large density perturbation for case (d). However there are some quantitative differences 
among three resolutions, cases (a) and (b) have milder steeps in the exponential growth of odd modes than cases (c) 
and (d). This might come from the diffusive effect of the large meshes of the computational grid. Hereafter, we 
will not discuss the precise growing feature, but rather discuss the mechanism of the odd azimuthal modes in the 
dynamical bar-instability. In this case, our choice of the computational resolution (case(b)) can explain the growth 
of odd azimuthal modes for the above purpose. The growth of the odd modes can be seen in both relativistic and 
present Newtonian results, and should be worthy of investigation. In order to study it by further numerical works, 
different approaches or much sophisticated treatments would be necessary. We here pay attention to a possible growth 
mechanism of the odd modes. 
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III. ONE-DIMENSIONAL FLOW COUPLED TO A SCALAR FIELD 

A. Mathematical Model 

It is clear that the most important nonlinearity in eqs.([T])-(j4]) is the advection term, (uV)v, when written in non- 
conservative form. Taking into account the nonlinearity and instability caused by an external force, we introduce a 
simplified model to examine the nonlinear evolution of the unstable modes. Our model is Burgers' equation for a flow 
velocity u(t, x) coupled to a scalar field <f>{t,x): 

d t u + ud x u — vd 2 u + X(j), (16) 

dl<f> + 2<j> = -u + l, (17) 

where v is a diffusion constant. It is well known that ea. (|16p with A = represents a shock model due to the nonlinear 
advection term ud x u. We assume v > for stability. The additional term \<p is a force mimicking gravity and may 
cause global instability, as shown below. The scalar field 4> satisfies a second-order partially differential equation (fl7|) 
which is like the Poisson equation with a source u. The additional term 2<j>, which may come from geometrical factors, 
such as curvature, in a realistic system, is introduced to adjust the m = 2 mode to the most unstable one. See the 
next subsection. We assume that the system of eqs. ([TBI and (fl~7|) is in non-dimensional form and that the spatial 
range is limited to < x < 27r. Periodic conditions are imposed on the functions u and <j> at x — 0, 27r. Therefore, x 
corresponds to azimuthal angle </? in the realistic system. It is easily checked that 1 — u = = Oisan exact solution 
of eqs. (HHJ) and (fTT)) . This solution corresponds to an axisymmetric solution for eqs.(UJ)-(|2]). We regard u = 1 and 
<f> = as the background state and consider linear stability and nonlinear growth from this uniform state. 
We solve eqs. (fTT))) and (fTT)) using Fourier series expansion for < x < 2tt: 

u = 1 + a m (t) cos(mx) + b m (t) sm(mx). (18) 

m— 1 

From eq. (fT7[) we have 

<p = 1 {a m (t) cos(mx) + b m (t) sin(mx)) . (19) 

m—l 

The number and growth rates of the unstable modes are easily controlled by changing the initial data and the 
parameters v and A. Our model is one-dimensional and is, therefore, easily solved for a wide range of parameters. 
We will show that this model's description of mode coupling, unstable growth and decay mimics the realistic system 
of eqs.([T])-([4]), very well. 

B. Linear Perturbation 

Assuming that |a m |, \b m \ <C 1, we linearize eq. ([T6")) giving 

^(a m ± ib m ) T im(a m ± ib m ) = T m (a m ± ib m ), (20) 

where 

T m = ~vm 2 + — ^r. (21) 
m z — 2 

The solution can be written as a m ± ib m cx exp(±imi + T rn t), and the stability of mode m is therefore determined 
by the sign of T m . That is, r m represents the growth rate (for T m > 0) or the decay rate (for r m < 0). It is easily 
seen that the diffusion term with v > is stabilizing, while the term Xcj) with A > is de-stabilizing. The growth 
rate depends on the magnitudes of A and u, whereas the number of unstable modes depends only on the ratio of two 
constants \jv. The m mode becomes unstable if \/v > m 2 (m 2 — 2). This model shows that as A increases, the short 
wavelength modes, i.e., those with large m, become unstable. For example, all modes are stable for \jv < 8. For 
8 < \ jv < 63 the only unstable mode is m = 2. For 63 < X/u < 224 the unstable modes are m — 2 and m = 3. 

However, the m = 1 mode is always stable. This choice of the model may correspond to some constraints in a 
realistic system, e.g., no motion of the center of mass. Furthermore, the most unstable mode is adjusted to m — 2 in 
this model. The correspondence to the realistic hydrodynamical system is clear with respect to wave number m, but 
we can not at present demonstrate that the system of eas. (|16p - (ll7|) is approximately derived from a realistic system. 
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C. Numerical Calculation for Nonlinear Growth 

In this section we use the Fourier series expansion (|18p . but do not assume that a m (t), b m (t) are small. From eq. (|16[) 
we have a coupled system of ordinary differential equations for a m (t) and b m (t) (m = 1,2, •• •). In the numerical 
calculations, the range of m is truncated to m max = 40. The parameter A is set to A = 1, so that the growth timescale 
is the dynamical one. The viscous timescale associated with v should be less than the dynamical one. The value v 
is quite small in realistic situations, but is not so small in the model calculation. We used v = 0.01 ~ 0.15 to save 
simulation time and control the number of unstable modes. For this choice of v, the viscous timescale ~ v m , 
given by eq. (|21[) . is larger than the dynamical one for small m modes, but not for large m modes, e.g, m > v~ x / 2 . 
The decay rate T m is, therefore, modified so as to suppress the very rapid decay of large m modes. The maximum 
decay rate is set to |T m | = 1. 

It is well known that some finite difference schemes for integrating the inviscid Burgers' equation, i.e., v = A = in 
eq.(|16p. become unstable at the shock[27|: overstable oscillations with high frequency are generated. Our numerical 
scheme is not a finite difference scheme but we tested it for the case v = A = and found the scheme to be stable. 
One drawback of the numerical method is the Gibbs phenomenon: an overshoot at the shock front. The oscillation 
cannot be removed by increasing the number m max . However, this peculiarity always appears at a discontinuity even 
for a static problem. For example, the discontinuity of a sawtooth form is poorly expressed by a Fourier sum[28l|. 
Furthermore, note that the Gibbs phenomenon occurs with other cigenfunction expansions. 

The time-evolution of the amplitudes log(C m ) = \og 10 (y/a 2 n + fe^J for some Fourier modes (m = 1, • • • , 6) is shown 
in FigurelH Figure[4ja) shows the results for X/u » 6.7 (A = 1, v — 0.15). The linearized perturbation theory predicts 
no growing modes for this system. The initial amplitudes are chosen as ai = = 10~ 2 and the others are zero. 
Higher m modes are always induced by the initial seeds, 02 or 03, but they all decay with time and their amplitudes 
are small compared with the m = 2 and m = 3 modes. That is, the m = 2 mode is the largest of the even modes and 
the m = 3 mode is the largest of the odd modes. 

Figure BJb) shows the results for a less viscous model with X/v = 20 (A = l,u = 0.05). There is one unstable 
mode, the m = 2 mode, predicted by linearized perturbation theory. The initial conditions are the same as those in 
Figure SJa). The m — 2 mode grows exponentially until t ~ 15, where the amplitude of the m — 2 mode reaches the 
nonlinear regime: 10~ 2 x exp(0.3 x 15) ~ 1. The growth rate until t ~ 15 agrees with T2 — 0.3 derived from the 
linearized theory. All other even modes, originating from the bilinear coupling term ud x u, also grow. The m = 6 mode 
is produced from the coupling between m = 2 and m = 4 and also from the quadric coupling of m = 3. Therefore, 
the amplitude of the m = 6 mode is not always smaller than that of m = 4. The growth of all even modes is slightly 
suppressed after the turning time t ~ 15. The turning time is also important for the odd modes. The odd modes 
decay for t < 15, but grow after that. Therefore, the nonlinearity of the amplitude of the m = 2 mode cannot be 
ignored even for the odd modes. The turning time corresponds to shock formation as will be discussed later. In order 
to examine the effect of the nonlinearity of the m — 2 mode on the growth of all other modes, we set the initial 
amplitude of 02 to 10~ 4 . The other initial conditions were kept the same as those used to produce Figure 0Jb). The 
time-evolution shown in Figure \M,c) has the same general features as Figure BJb) , but the turning time, when the 
odd modes switch from decay to growth, is shifted to t ~ 30. This is because of the small initial amplitude of the 
?7i = 2 mode: 10 -4 x exp(0.3 x 30) ~ 1. For eqs. tfTrj)) and (|17[) all odd modes are always zero, if they are exactly zero 
initially. When there is at least one odd mode with a finite amplitude, the nonlinearity of the m = 2 mode enhances 
all odd modes. 

The nonlinear evolution for X/v = 100 (A = 1, v = 0.01) is shown in FigureBJd). The initial conditions are the same 
as those used in Figure Bib). In this model the m = 2 and m = 3 modes are unstable with growth rates of T2 = 0.46 
and T3 = 0.053 from the linearized theory. Overall the features are the same as in Figure|U[b), except for the timescale. 
The turning time due to the m = 2 mode is shorter in this model: t ~ 10 since 10~ 2 x exp(0.46 x 10) ~ 1. The m = 3 
mode is initially unstable but does not grow significantly during the early phase t < 10. The typical growth timescale 
is very long, T^" 1 ~ 20, so that mode coupling becomes much more important at early times. However, the unstable 
m = 3 mode maintains the amplitudes of other odd modes at higher levels through mode coupling before the turning 
time t ~ 10. 



D. Comparison with 3D Simulations 

By comparing the mathematical model with 3D numerical results, the following features become clear. The odd 
modes grow only after nonlinear saturation of unstable mode. This fact can easily be seen in Figure [21 in which the 
initial amplitude of odd modes is 10~ 3 and the nonlinearity becomes important in short timescale. We further discuss 
the growth feature in the models with small initial amplitudes. In the mathematical model, the odd mode growth 
starts approximately from the time t s + 2r, where t s and t are the saturation time and growth time of the unstable 
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(a) (b) 




FIG. 4: The time-evolution of the amplitudes of the lowest six Fourier modes (m = 1, • • • 6). Panel (a) is result for v = 0.15, 
(b) and (c) for v = 0.05, and (d) for v — 0.01. The initial amplitudes are ai = az = 10~ 2 for (a), (b) and (d), while 
a,2 = 10 _4 ,a3 = 10 -2 for (c). The behavior of very small amplitude, say, log(C m ) approximately less than —10 in panel (a), 
mostly comes from numerical truncation errors, and is therefore unimportant. 



mode. The starting point of the odd mode growth is not easily determined in the actual calculations, since the growth 
curve is not so sharp. The similar relation is however realized. The amplitudes of odd modes exceed, say, ~5x 10~ 3 
at the time ~ t s + 5r in the models I-III in Figure [TJ This property is almost independent of initial perturbations 
with 5 < 10 -3 and unstable models. The start time t g of odd mode growth is roughly given by t g = t s + err, where 
a = 0(1). The number a depends on dynamical degrees of freedom of the system. The growth of the amplitude 
a m simply depends on the bilinear coupling Y^o>kO"m-k or J2 a kb m ~k in the mathematical model. There is a similar 
bilinear coupling in 3D simulation, but coupling is more complicated. The azimuthal Fourier component of the density 
couples with three components of velocity, which couple with those of energy, gravity and density. Moreover, they 
are functions of x,y,z. The dynamical degrees are so large in the 3D system, that the interval (t g — t s )/r becomes 
longer. The relation between the nonlinear saturation time and growth time of odd modes qualitatively holds in 
both systems, despite of the increased dynamical degrees of freedom. In the mathematical model, it is possible to 
eliminate the initial perturbations of odd modes. In such a clean case, the odd modes can not appear. In the actual 
numerical simulations, some random noises, whose amplitudes are expected as ~ 10~ 5 -10~ 4 , should be involved. The 
noise reduction would be possible in future simulation on the much large scale computer, but is not necessary in the 
realistic applications since such an ideal initial condition is rare. Thus odd modes would appear in general. 
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FIG. 5: Snapshots of the function u(t,x) are shown as a function of x' = x/2n for t = 0, Air, 87r, 167r. The attached labels 
denote the time 

E. Evolution of Shape 

The similarity can be seen in the time evolution of the Fourier components both in mathematical and 3D numerical 
models as shown in previous subsection. The time evolution of the shape u(t, x) is shown in Figure[5] The parameters 
and initial conditions are the same as used to produce Figure H^b). The snapshots are given for times t — 0,47r, 87r 
and 167r. The choice comes from removing the propagation effect, because the initial velocity is u ss 1. The m = 2 
mode initially grows and the shape is enhanced before the turning time t ~ 15. The curve at t = An clearly shows 
symmetric features due to the m = 2 mode. That is, the shape is the symmetry under translations x — > x + tt, a 
"7r-symmetry" . The nonlinearity causes a shock as in the original Burgers' equation. After shock formation, the Gibbs 
phenomenon associated with Fourier series is seen at t — 8tt, 16ir. The overshoot is a numerical artifact and such 
behavior always appears when a function having a sharp discontinuity is expressed as a Fourier series [28j|. Neglecting 
the Gibbs phenomenon, the symmetry due to the m = 2 mode can still be seen in the shape at t = 8tt, whereas it 
is partially broken at t — 167r. The time t = 1&tt in the mathematical model is much longer than that of non-linear 
saturation and that of growth of odd modes. Therefore, there is no counterpart in 3D numerical simulations in Figures 
1 and 2. The mathematical model suggests that a "7r-symmetry" (i.e., symmetry under a 180° rotation around the 
z-axis) in the shape is broken in a longer timescale. 



IV. DISCUSSION 



We have considered the nonlinear evolution of the bar-type instability in a differentially rotating star with significant 
rotational energy. The previous numerical results of general relativistic simulation 14) suggest that the growth is likely 
to come from mode coupling. In order to obtain further evidence of mode coupling, we have developed a simulation 
of three-dimensional hydrodynamics in Newtonian gravity and a simple mathematical model. Our mathematically 
simplified model provides a concrete example showing the importance of mode coupling. The amplitudes of odd 
modes increase without unstable odd modes being present in the axially symmetric state; instead, they are enhanced 
by the bar instability with m = 2. We also confirmed that this physical picture is consistent with the results from a 
three-dimensional hydrodynamics simulation. Generally, the odd modes grow only after the bar instability reaches the 
nonlinear regime. The timescales of the mode coupling and the growth of unstable modes may depend on the rotation 
law and the strength of the initial instabilities. It is very rare that the initial perturbations in the hydrodynamics 
simulation should consist of purely even or odd modes only. Therefore, the unstable bar mode enhances the amplitudes 
of the all other modes at late times, no matter whether they are even or odd. 

A similar mode coupling can be seen in numerical simulations for the one-armed spiral instability [l7j and the 
elliptical instability[l9j of rotating stars in Newtonian gravity. The initial models and the growth mechanism are 
different, but the turbulent-like behavior appears in diagnostics of the azimuthal Fourier components at late times of 
nonlinear growth [IH [l9j. The behavior is also important for the nonlinear saturation of the unstable mode. Further 
study is necessary to explore the origin of the similarity seen in the development of different unstable modes. It is 
reasonable to assume that the nonlinearity in hydrodynamics is the source of this similarity. What is the effect of 
general relativity? A number of nonlinearities occur in general relativity which may affect the growth of the unstable 
bar mode. Although the time-evolution in full relativistic calculations is very similar to that in Newtonian gravity [l4j. 
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it will be very interesting to explore further whether or not a full relativistic simulation produces a nonlinearity different 
from the one presented by the simple model. 
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